function G = slvecgraph(X, X2, edges, f, varargin)
%SLVECGRAPH Creates a graph among vectors
%
% [ Syntax ]
%   - G = slvecgraph(X, [], edges, f, ...)
%   - G = slvecgraph(X, X2, edges, f, ...)
%
% [ Arguments ]
%   - X:        the (source) sample matrix
%   - X2:       the target sample matrix
%   - edges:    the matrix representing the edge set (ne x 2)
%   - f:        the function to compute edge values
%
% [ Description ]
%   - G = slvecgraph(X, [], edges, f, ...) creates a graph among
%     samples by using a specified function handle to compute edge
%     values.
%
%     Suppose there are n samples in X, then X should be a d x n 
%     sample matrix, with each column representing a sample vector.
%     This function constructs a graph with n nodes, and each node
%     corresponds to a sample.
%
%     If there are ne rows in edges, then edges(i, :) represents the
%     i-th edge. G is the adjacency matrix of the graph. If there is 
%     an edge connecting the i-th node and j-th node, then G(i, j) equals
%     the corresponding edge value, otherwise G(i, j) is zero.
%
%     The edge values are computed using the function handle f, which
%     is invoked with the following syntax
%       $ evals = f(V1, V2) $
%     V1 and V2 are two subset of sample vectors, and they contain the same
%     number of samples, say n. Then evals should be a row vector of size
%     1 x n, in which evals(k) is the value corresponding to the edge
%     connecting the sample vectors V1(:,k) and V2(:,k).
%
%   - G = slvecgraph(X, X2, edges, f, ...) creates a bigraph, of which the
%     source nodes correspond to the vectors in X, while the target nodes
%     correspond to the vectors in X2. 
%
%     If X and X2 respectively have n and n2 samples, then G will be a
%     matrix of size n x n2.
%
%     The following options can be specified for graph construction.
%     \{:
%         - output:     the output form of the adjacency matrix G
%                       (default = 'sparse')
%                       \{:
%                           - full:     G is created as a full matrix
%                           - sparse:   G is created as a sparse matrix
%                       \:}
%
%         - vcahce:     the maximum size (in term of the number of scalars)
%                       of vectors stored in value-computation-cache.
%                       (default = 1e7).
%     \:}
%
% [ History ]
%   - Created by Dahua Lin, on Jul 7, 2007
%

%% parse and verify input arguments

error(nargchk(4, inf, nargin));

assert(isnumeric(X) && ndims(X) == 2, 'sltoolbox:slvecgraph:invalidarg', ...
    'The X should be a 2D sample matrix');
if isempty(X2)
    X2 = X;
else
    assert(isnumeric(X2) && ndims(X2) == 2, 'sltoolbox:slvecgraph:invalidarg', ...
        'The X should be a 2D sample matrix');
end

n = size(X, 2);
n2 = size(X2, 2);
pairlen = size(X, 1) + size(X2, 1);

assert(isnumeric(edges) && ndims(edges) == 2 && size(edges, 2) == 2, ...
    'sltoolbox:slvecgraph:invalidarg', 'The edges should be a numeric array with two columns');
ne = size(edges, 1);

assert(isa(f, 'function_handle'), 'sltoolbox:slvecgraph:invalidarg', ...
    'f should be a function handle.');

opts = struct( ...
    'output', 'sparse', ...
    'vcache', 1e7);
opts = setopts(opts, varargin{:});

assert(ischar(opts.output), 'sltoolbox:slvecgraph:invalidopt', ...
    'The option output should be a string');
opt_output_idx = find(strcmp(opts.output, {'sparse', 'full'}));
assert(~isempty(opt_output_idx), 'sltoolbox:slvecgraph:invalidopt', ...
    'The option output should be either ''sparse'' or ''full''.');
is_sparse = (opt_output_idx == 1);

assert(isnumeric(opts.vcache) && isscalar(opts.vcache) && opts.vcache > 0, ...
    'sltoolbox:slvecgraph:invalidopt', 'The option vcache should be a positive scalar value.');
assert(opts.vcache > pairlen, 'sltoolbox:slvecgraph:invalidopt', ...
    'The option vcache should be large enough to hold at least a pair of sample vectors.');

%% main 

% allocate graph matrix

if is_sparse
    G = spalloc(n, n2, ne);
else
    G = zeros(n, n2);
end

% compute and set edge values
s_idx = edges(:, 1);
t_idx = edges(:, 2);

einds = sub2ind([n, n2], s_idx, t_idx);

if pairlen * ne < opts.vcache % single batch
    evals = f(X(:, s_idx), X2(:, t_idx));
    
else    
    bn = floor(opts.vcache / pairlen);
    sp = 1:bn:ne;
    ep = [sp(2:end), ne];
    nb = length(sp);
    
    evals = zeros(ne, 1);
    for i = 1 : nb
        csp = sp(i);
        cep = ep(i);
        
        evals(csp:cep) = f( X(:, s_idx(csp:cep)), X2(:, t_idx(csp:cep)) );
    end
    
end
  
G(einds) = evals;

